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1.0  SUMMARY 


This  final  report  describes  the  effort  to  develop  a  highly  accurate  and  efficient  simulation  tool  for 
the  multi-physics  and  multiscale  modeling  and  simulation  of  the  electromagnetic-plasma 
interaction  and  the  high-power  microwave  breakdown  in  air.  Under  the  high  pressure  and  high 
frequency  condition  of  the  high-power  air  breakdown,  the  physical  phenomenon  is  described 
using  a  nonlinearly  coupled  full-wave  Maxwell  and  fluid  plasma  system.  This  nonlinear  system 
is  then  solved  by  a  numerical  method  that  is  uniformly  high  order  in  both  space  and  time. 
Accurate  and  efficient  numerical  simulation  methods  and  techniques  are  introduced,  developed, 
and  validated.  Several  numerical  examples  are  given  to  demonstrate  the  capability  of  the 
numerical  method  in  the  modeling  and  simulation  of  the  high-power  microwave  air  breakdown 
problems,  where  the  underlining  physical  process  can  be  well  understood  and  interpreted  from 
the  simulation  results. 

2.0  INTRODUCTION 

High-power  microwave  (HPM)  devices  and  systems  have  very  important  civilian  and  military 
applications.  To  design  better  HPM  devices  that  generate  higher  electromagnetic  power  and 
longer  pulse  width,  extensive  research  effort  has  been  devoted  to  the  development  of  microwave 
sources  [1][2],  the  design  of  output  windows  [3][4],  and  the  optimization  of  advanced  cathodes. 
However,  as  the  power  density  goes  higher  and  higher,  the  HPM  devices  and  systems  are  more 
and  more  vulnerable  to  the  HPM  breakdown,  including  the  air  breakdown  during  the  generation 
and  transmission  of  HPM  and  the  dielectric  window  breakdown  when  the  HPM  is  radiating 
through  the  dielectric  window.  When  breakdown  happens,  the  transmission  capability  of  a 
microwave  device  can  be  severely  limited. 

The  HPM  breakdown  can  occur  inside  the  HPM  device  and  on  the  output  window  due  to  surface 
flashover  phenomena  [5]-[7]  and  multipactor  [8][9],  During  the  HPM  generation  and  radiation,  the 
electrons  in  the  HPM  device  are  accelerated  by  the  high-intensity  electromagnetic  fields.  The  motion 
of  electrons  produces  impact  on  other  gas  particles  or  metal  and  dielectric  surfaces,  which  can  release 
secondary  electrons  into  the  device.  These  electrons  are  then  accelerated  by  the  electromagnetic 
fields  and  also  impact  other  particles  or  surfaces.  Meanwhile,  they  also  generate  secondary 
electromagnetic  fields  and  enhance  the  original  fields.  If  the  power  density  is  sufficiently  high,  the 
free  electrons  can  sustain  themselves  through  the  electromagnetic  fields  they  generate,  and  therefore, 
result  in  an  exponential  multiplication  of  the  electrons,  which  is  known  as  breakdown.  This  may  lead 
to  the  malfunction  or  even  damage  of  the  HPM  devices  and  systems.  To  achieve  a  better 
understanding  of  the  breakdown  mechanism  and  process  and  to  improve  the  design  of  microwave 
devices,  research  investigation  has  been  conducted,  both  experimentally  [3] [4]  and  computationally 
[10][1 1],  By  comparing  the  numerical  results  obtained  from  numerical  simulations  with  experimental 
observations,  breakdown  theories  can  be  validated  and  the  underlying  physics  can  be  better 
understood.  In  order  to  simulate  such  a  multi-physics  process  involving  field-particle  interaction  and 
strong  nonlinearity  with  a  high  accuracy  and  reliability,  both  the  electromagnetic  physics  and  the 
plasma  physics  have  to  be  described,  modeled,  and  solved  accurately. 

Governed  by  Maxwell’s  equations,  the  electromagnetic  fields  can  be  solved  with  various 
numerical  methods  developed  by  the  computational  electromagnetics  community.  Among  all  the 
methods,  the  finite-difference  time-domain  (FDTD)  method  is  most  widely  used  because  of  its 
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simplicity  and  high  efficiency.  Also  because  of  its  scalar  representation  of  each  field  component 
on  structured  mesh  grids,  the  continuity  of  the  electromagnetic  fields  can  be  preserved,  which  is 
especially  important  when  the  electromagnetic  fields  are  coupled  to  the  particle  motion. 
However,  the  stair-case  approximation  of  the  solution  domain,  the  finite-difference 
approximation  of  the  fields  and  their  derivatives,  and  the  explicit  leap-frog  time-marching 
scheme  used  in  the  FDTD  method  result  in  a  low-order  accuracy,  which  requires  an  extremely 
dense  mesh  grid  and  an  extremely  tiny  time  step  size  in  a  simulation  to  achieve  a  desired 
accuracy.  On  the  contrary,  another  widely  used  method,  the  finite-element  time-domain  (FETD) 
method  [12],  which  employs  an  unstructured  mesh  for  the  geometric  discretization,  high-order 
vector  basis  functions  for  the  field  expansion,  and  implicit  time-marching  scheme,  is  able  to 
achieve  high-order  accuracy  and  unconditional  stability.  Unfortunately,  the  implicit  scheme 
results  in  a  globally  coupled  system  that  has  to  be  solved  at  every  time  step,  which  leads  to  a 
reduced  efficiency  compared  to  the  FDTD  method.  More  importantly,  because  of  the  application 
of  the  vector  basis  functions,  only  the  field  components  that  are  tangential  to  the  elemental 
interfaces  are  continuous,  and  those  that  are  normal  to  the  interfaces  are  much  less  accurate. 
When  coupled  with  the  particle  solver,  the  discontinuous  normal  components  of  the  fields  will 
impose  an  unpredictable  amount  of  force  to  the  particles  and  result  in  a  spurious  solution  or  even 
a  numerical  breakdown  of  the  simulation.  To  overcome  the  abovementioned  issues  of  both 
methods,  in  this  work  the  nodal  discontinuous  Galerkin  time-domain  (DGTD)  method  [13]-[19] 
is  employed  to  solve  for  the  electromagnetic  fields.  Different  from  the  vector  DGTD  [20]  [21] 
traditionally  used  in  the  computational  electromagnetics  community,  the  nodal  DGTD  method 
employs  nodal  interpolatory  basis  functions  to  represent  the  physical  quantities,  which 
guarantees  the  field  continuity  in  all  three  spatial  directions.  Using  an  unstructured  mesh,  high- 
order  nodal  basis  functions,  and  an  explicit  high  order  time-marching  scheme,  the  nodal  DGTD 
method  is  able  to  achieve  high-order  accuracy  in  the  geometric  discretization  of  the  solution 
domain  and  in  the  spatial  and  the  temporal  discretization  of  the  field  solution.  In  addition,  the 
nodal  DGTD  method  also  enjoys  the  following  advantages:  (1)  high  parallel  efficiency  on  a 
massively  parallel  platform;  and  (2)  high  flexibility  and  efficiency  by  using  dynamic  /^-adaption 
and  local  time-stepping  scheme. 

The  other  half  of  the  physical  phenomenon,  the  plasma  physics  [22]  [23],  can  be  described  by  the 
Vlasov  and  the  Boltzmann  equations  in  the  collision-less  and  the  collisional  cases,  respectively. 
Because  of  the  high  dimensionality  of  the  Vlasov  and  the  Boltzmann  equations  and  the 
consequent  difficulty  in  their  direct  solution,  various  simplified  models  are  usually  used,  where 
in  the  low  to  moderate  pressure  and  density  case,  the  well-known  particle-in-cell  (PIC)  method 
has  been  widely  adopted  in  the  simulation  of  high-power  microwave  devices,  and  in  the  high 
pressure  and  density  case,  fluid  models  work  adequately  to  describe  the  plasma  behavior.  In  the 
case  of  air  breakdown  under  atmosphere  condition,  fluid  models  can  be  adopted.  By  taking  the 
first  several  moments  of  the  Vlasov  or  the  Boltzmann  equation,  a  hierarchy  of  models  can  be 
obtained,  from  the  simplest  two-moments  drift-diffusion  model  [24],  to  the  three-moments 
hydrodynamic  model  [25],  to  the  four-moments  energy  transport  model  [26] [27],  to  the  six- 
moments  model  [28],  and  so  on.  With  the  inclusion  of  higher  order  moments,  more  physical 
mechanisms  are  accounted  for  by  the  more  sophisticated  model.  These  models  govern  plasma 
quantities  such  as  the  plasma  density,  velocity,  temperature,  and  energy.  Since  the  related  plasma 
parameters,  such  as  ionization  frequency,  collision  frequency,  attachment  frequency,  and 
electron  energy  loss  frequency,  become  highly  nonlinear  when  the  electromagnetic  field  intensity 
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is  sufficiently  high,  the  resultant  coupled  system  also  becomes  highly  nonlinear.  Such  a  strong 
nonlinearity  experienced  in  HPM  devices  inhibits  the  current  computational  methods  from 
simulating  the  breakdown  process  accurately  and  efficiently.  As  a  result,  the  major  techniques 
used  to  increase  the  breakdown  threshold  in  HPM  devices  and  systems  are  still  empirical. 

In  this  work,  the  nodal  DGTD  is  used  to  simulate  both  the  electromagnetic  physics  and  the 
plasma  physics.  In  Section  3.0,  the  physical  model  and  its  challenge  to  the  numerical  methods 
will  first  be  introduced,  and  the  proposed  numerical  scheme  and  techniques  will  be  described  in 
detail.  Several  numerical  examples  will  be  presented  in  Section  0,  and  conclusion  will  be  drawn 
in  Section  0. 


3.0  METHODS,  ASSUMPTIONS,  PROCEDURES 

3.1  Physical  Model  and  Its  Challenges 

3.1.1  Plasma  Fluid  Model 


The  plasma  model  adopted  in  this  work  is  a  simple  drift-diffusion  model  that  was  originally 
developed  in  [29]-[33].  Starting  from  the  Boltzmann  equation  (general  form) 


df  F  (df\ 

-UvV/  +  --Vv/=  -M 

Ot  TTL  Vtft/coll 


(i) 


where  /  =  f(r,v,  t)  stands  for  the  seven-dimensional  probability  distribution  function  (PDF)  of 
a  single  particle  species,  with  the  seven-dimensional  phase  space  spanned  by  (r,  v,  t),  F  stands 
for  the  force  presented  to  the  particle  with  a  mass  m,  the  term  on  the  right-hand-side  (RHS) 
stands  for  the  collision  mechanism.  In  the  fluid  model,  the  macroscopic  concepts  are  defined  by 
integrating  the  PDF  over  the  velocity  dimensions,  which  yields  the  definition  of  the  particle 
density 


the  mean  velocity 


f(r,v,  t)  dv 


u  = 


vf(r,v,t )  dv 


and  the  mean  energy 


m(  .  m  f 

£  =  y  >  =  y^  I  v  f(r,v,t )  dv 
of  the  particle  species  under  consideration. 


(2) 


(3) 


(4) 
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The  basic  conservation  laws  of  the  plasma  physics  can  be  derived  by  taking  the  first  three 
moments  of  the  Boltzmann  equation,  which  results  in  the  mass  conservation  equation  (aka  the 
particle  continuity  equation) 


dn 

—  +  V  ■  (nu)  =  5 
dt 


with  5  being  the  source  term,  the  momentum  conservation  equation 

dnu  1  -  F 

— — I-  V  ■  (nuu)  = - V  ■  P  +  n  — l-i? 

dt  mm 

and  the  energy  conservation  equation 

dn£ 


dt 


+  V  ■  (nuS  +  P  ■  u  +  Q)  =  nu  ■  F  +  S£ 


(5) 


(6) 


(7) 


respectively. 

In  the  case  of  high-power  microwave  breakdown  in  air,  the  following  reasonable  assumptions 
can  be  made.  (1)  The  pressure  tensor  is  P  is  diagonal  and  isotropic,  hence,  is  a  scalar.  (2)  Since 
the  discharge  takes  place  at  high  pressure,  the  high  collision  condition  can  be  assumed.  (3)  Local 
equilibrium  between  electric  acceleration  (energy  gain  from  the  electric  field)  and  collisional 
momentum  and  energy  loss  is  assumed,  so  that  the  ionization  frequencies  depend  only  on  the 
local  electric  field.  This  is  known  as  the  local  field  approximation.  Under  these  assumptions,  one 
can  end  up  with  the  drift-diffusion  equation 

dn 

—  +  V  ■  (±nnE  -  DVn)  =  S  (8) 

dt 

where  /r  is  the  mobility  coefficient  and  D  is  the  diffusion  coefficient.  Since  the  breakdown  under 
consideration  takes  place  at  microwave  frequencies,  the  drift  term  zeros  out  over  a  full  period  of 
the  electromagnetic  wave,  and  therefore,  can  be  ignored.  As  a  result,  the  diffusion  mechanism 
dominates  over  the  drift  mechanism  in  the  microwave  frequency  range.  Under  the  quasi-neutral 
assumption  of  the  plasma,  and  taking  into  account  of  the  ionization,  attachment,  and 
recombination  effects,  the  drift-diffusion  equation  finally  becomes 

dn 

—  =  V  ■  (DgffVn)  +  (u;  -  ua)n  -  rein2  (9) 

where  Deff  is  an  effective  diffusion  coefficient,  u;  and  va  are  the  ionization  and  the  attachment 
frequencies,  respectively,  and  rei  is  the  electron-ion  recombination  coefficient. 

In  the  above  equation,  the  effective  diffusion  coefficient  Deff  is  a  combination  of  the  free  electron 
diffusion  coefficient  De  and  the  plasma  bulk  ambipolar  diffusion  coefficient  Da 
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(10) 


Dcff  — 


aDe  +  Da 
a  +  1 


where  a  is  a  weighting  parameter  given  by 


a  =  UjTM  =  AB/L2 


(11) 


with  tm  =  €0/[en(/2e  +  /r;)]  standing  for  the  local  dielectric  (Maxwell)  relaxation  time,  e0  being 
the  vacuum  permittivity,  and  /te  =  e/mevm  and  =  /ic/200  being  the  electron  and  ion 
mobility,  respectively.  The  local  Debye  length  AD  and  the  characteristic  diffusion  length  L  are 
given  by  (kB  is  the  Boltzmann  constant,  and  Te  is  the  electron  temperature) 

Ad  =  yfe^Tjehi  (12) 


and 


L  =  ||Vn/n||  1  =  (13) 

respectively.  Around  the  edge  of  the  plasma  bulk,  since  the  electron  density  drops  quickly  from  a 
high  level  to  almost  zero,  the  electrons  in  this  region  diffuse  with  the  electron  free  diffusion 
coefficient 


Dc  =  kc 


e 


(14) 


In  the  plasma  bulk,  because  of  the  presence  of  both  electrons  and  ions,  the  much  slower 
ambipolar  diffusion  is  the  dominant  mechanism.  Since  the  ions  are  much  heavier  than  the 
electrons,  the  ion  diffusion  coefficient  is  very  small  and  is  usually  ignored 

kBT{ 

^  —  *  0.  (15) 

*7i 


The  ambipolar  diffusion  coefficient,  as  a  combination  of  the  electron  and  the  ion  diffusion 
coefficients,  is  hence  given  by 


D 


a 


^De  +  iieDt  ^  kBTe 

- ~  —  D„  =  n, - 

ki  +  kc  Me  e 


(16) 


since  /re  »  /ij  and  Dc  »  Dt.  To  calculate  both  the  electron  and  the  ambipolar  diffusion 
coefficients,  the  electron  temperature  needs  to  be  obtained,  which  is  usually  expressed  by  an 
empirical  relation  between  the  electron  temperature  Te  and  the  effective  field  £cff  (will  be  defined 
in  Section  3.2) 
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(17) 


l 


_ 

2.1  x  10"3  — | 

(  FcffY 

91.0 +  —) 

e 

P 

V  P  J\ 

Other  important  physical  mechanisms  include  ionization,  attachment,  and  recombination.  The 
ionization  and  attachment  frequencies  are  usually  combined  as  the  effective  ionization  frequency 
ueff  which  can  be  expressed  in  terms  of  the  electron  drift  velocity  vd 


veff  =  Vi  -  va  =  avd  (18) 

vd  =  peE  (19) 


where  a  is  the  ionization  coefficient  which  characterizes  the  number  of  ionization  events  that  an 
electron  undergoes  per  unit  length  along  the  field  and  is  given  by  the  following  empirical 
expression 


a  =  A0p[ex p-Bo(p/£eff-p/£c)  _  i]  (20) 

if  £eff/p  <  50  (V/cm.torr),  and 

a  =  Apexp~Bp/,Ee{i  (21) 

if  Eeff/p  >  50  (V/cm.torr).  In  the  above  expressions,  A0  =  0.005  (cm_1torr_1),  B0  =  200 
(V/cm.torr),  and  Ec/p  =  31.25  (V/cm.torr),  where  Ec  stands  for  the  critical  field  intensity 
beyond  which  the  ionization  will  dominant  over  attachment.  In  air,  the  coefficients  A  =  8.805 
(cm_1torr_1),  B  =  258.45  (V/cm.torr)  when  50  <  Eeff/p  <  200  (V/cm.torr),  and  A  =  15 
(cm_1torr_1),  B  =  265  (V/cm.torr)  when  100  <  Eeij/p  <  800  (V/cm.torr). 

To  complete  the  description  of  the  plasma,  another  equation  that  governs  the  particle  velocity  is 
needed.  From  the  momentum  conservation  equation,  by  assuming  that  the  distance  travelled  over 
one  field  period  is  small  compared  to  the  length  scale  of  the  field  and  pressure  variation,  one  can 
obtain  the  kinetic  equation 


du 

dt 


(22) 


where  q  is  the  charge  carried  by  a  particle,  and  um  is  the  momentum  transfer  electron-neutral 
collision  frequency,  which  is  often  assumed  to  have  a  simple  relation  with  the  pressure,  for 
example 


—  =  5.3  X  109  s  1torr  1  (in  air). 
V 

3.1.2  Coupled  Electromagnetic-Plasma  System 


(23) 


The  interaction  between  the  electromagnetic  waves  and  the  plasma  fluid  is  governed  by  the 
coupled  electromagnetic-plasma  system,  where  the  electromagnetic  physics  is  described  and 
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governed  by  Maxwell’s  equations,  while  the  plasma  physics  is  described  and  governed  by  the 
electron  diffusion  equation  with  the  electron-neutral  collision,  ionization,  attachment,  and 
recombination  mechanisms  taken  into  consideration.  With  the  conduction  current  expressed  as 
/c  =  enu,  the  coupled  system  is  given  by 

1 


H  =  --V  x  E 

(24) 

p 

E  =  -(VxH-Jc ) 

(25) 

h=  »%/c+  m  («  +  «“) 

(26) 

h  =  V  ■  (DeffVn)  +  (i)[  -  ua)n  -  rein2. 

(27) 

Here,  the  scattered  field  formulation  for  Maxwell’s  equations  is  used  in  order  to  achieve  a  better 
numerical  accuracy  in  the  calculation  of  the  secondary  field  radiated  by  the  plasma  fluid,  which 
means  the  E  and  H  in  the  formulation  stand  for  the  secondary  (or  scattered)  field,  with  the 
superscript  “sea”  omitted  for  simplicity. 

3.1.3  Multi-Physics  and  Multiscale  Challenges 

Physically  speaking,  the  system  described  by  the  four  governing  equations  is  multiscale  in 
nature.  The  multiscale  nature  of  this  system  is  not  only  in  space  but  also  in  time.  To  understand 
such  a  multiscale  problem,  the  characteristics  of  the  electromagnetic  and  plasma  physics  are 
compared  by  taking  into  account  the  representative  values  of  the  physics  involved  in  a  typical 
application  scenario. 

For  the  spatial  domain,  the  characteristic  length  of  the  electromagnetic  physics  is  the  wavelength, 
which  is  on  the  order  of  millimeter,  for  example,  3  mm  at  100  GHz.  The  characteristic  length  of 
the  plasma  fluid  is  the  Debye  length  AD,  which  can  be  as  small  as  sub-micrometers.  For  example, 
for  the  electron  plasma  with  an  electron  temperature  of  2.0  eY,  and  an  electron  density  on  the 
order  of  1021/m3,  the  Debye  length  is  on  the  order  of  0.1  pn n. 

For  the  temporal  domain,  the  characteristic  time  of  the  electromagnetic  physics  is  the  period, 
which  is  on  the  order  of  picosecond,  for  example,  10  ps  at  100  GHz.  The  characteristic  time  of 
the  plasma  fluid  is  determined  by  the  diffusion  and  ionization  mechanisms,  which  is  on  the  order 
of  nanoseconds.  For  example,  under  the  same  assumption  for  the  electron  temperature,  and 
assume  the  magnitude  of  the  incident  electric  field  E  =  5  MY/m,  and  the  atmosphere  condition 
p  =  760  torr,  for  the  plasma  front  (considered  as  free  electrons)  to  propagate  a  distance  of  the 
characteristic  diffusion  length  L  =  ||Vn/n||-1,  at  the  estimated  propagation  velocity  v  = 
2 yfDgVi,  the  characteristic  time  is  estimated  to  be  on  the  order  of  1  ns. 

The  comparisons  made  here  are  summarized  in  Table  1.  Clearly,  for  the  spatial  domain,  the 
plasma  fluid  varies  faster,  but  for  the  temporal  domain,  the  electromagnetic  wave  varies  faster. 
To  design  a  good  numerical  algorithm  that  can  solve  the  strongly  coupled  nonlinear 
electromagnetic-plasma  system,  both  multiscale  characteristics  in  space  and  time  need  to  be 
considered  and  well  addressed. 
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Table  1.  Characteristics  of  the  Multiscale  Physics 


Electromagnetics 

Plasma  Fluid 

Ratio 

Fast  Physics 

Operating 

condition 

/  =  100  GHz 

p  =  760  torr 
kBTe  =  2.0  eV 
n  ~  1021/m3 

— 

— 

Spatial 

characteristic 

1  =  3  mm 

Ad  =  0.1  pm 

104: 1 

Plasma 

Temporal 

characteristic 

T  =  10  ps 

T  =  Ins 

1: 102 

Electromagnetics 

From  the  above  analysis,  it  seems  that  two  different  grids  can  be  used  in  the  spatial  discretization 
of  the  two  physics,  where  one  coarser  grid  that  resolves  electromagnetic  wavelength  can  be  used 
for  Maxwell’s  equation,  and  one  denser  grid  that  resolves  the  Debye  length  can  be  used  for  the 
electron  diffusion  equation.  Unfortunately,  it  is  not  the  case  here.  The  reason  is  that  the  spatially 
fast  varying  electron  density  distribution  and  its  resulting  electron  plasma  current  will  produce  a 
set  of  secondary  electric  and  magnetic  fields  that  are  varying  on  the  similar  scale  of  the  electron 
density  distribution,  which  is  much  faster  than  the  spatial  variation  of  the  incident  fields.  This 
can  be  observed  very  clearly  from  the  figures  in  Figure  1,  which  depicts  a  typical  application 
scenario  where  a  110-GHz  high-power  microwave  incidents  upon  and  interacts  with  a  plasma 
bulk  which  has  a  Gaussian  density  distribution  with  a  maximum  value  of  1015/m3  and  a 
standard  deviation  of  50  /on.  Notice  that  in  Figure  1(a),  the  wavelength  is  well  resolved  by  the 
discretization  grid  depicted  by  white  lines.  However,  with  the  same  grid,  the  zoom-in  plots 
shown  in  Figure  l(b)-(d)  reveal  the  fact  that  the  electron  density,  the  secondary  electric  and 
magnetic  fields  are  varying  on  the  same  spatial  scale,  and  hence  need  to  be  resolved  by  the  same 
resolution. 
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(C)  (d) 

Figure  1.  Spatial  Variation  of  the  Incident  Electric  Field  (a),  Initial  Electron  Density 
Distribution  (b),  Secondary  Electric  (c)  and  Magnetic  (d)  Fields  Produced  by  the  Electron 
Plasma  Current.  Only  %  of  the  Solution  Domain  (A  Spherical  Region  Filled  With  Air)  Is 

Shown  Because  of  Symmetry 


3.2  Numerical  Scheme 

In  this  section,  the  general  numerical  scheme  is  described,  with  all  the  mathematical  details  left 
for  the  next  section. 

With  the  physical  nature  and  challenges  of  the  problem  well  understood,  a  numerical  scheme  is 
designed  for  this  problem.  To  begin,  one  single  set  of  spatial  grid  is  used  for  both 
electromagnetic  and  plasma  physics.  Since  the  grid  cannot  be  extremely  small  to  resolve  the 
Debye  length,  higher  order  basis  functions  will  be  used  in  the  numerical  scheme  to  capture  the 
fast  variation  of  the  physical  quantities  within  each  of  the  grid  element.  The  use  of  the  higher 
order  basis  functions  can  achieve  higher  order  spatial  accuracy  without  a  dramatic  increase  of  the 
total  number  of  degrees  of  freedom  (DoFs),  which  sets  the  numerical  scheme  used  in  this 
research  apart  from  the  FDTD  method  used  in  the  past. 

After  the  spatial  discretization  issue  is  addressed,  the  four  partial  differential  equations  (PDEs) 
that  govern  the  coupled  electromagnetic-plasma  system  are  converted  into  four  ordinary 
differential  equations  (ODEs)  which  can  be  integrated  simultaneously  in  time  to  obtain  the  time- 
domain  solution.  However,  because  the  aforementioned  temporal  multiscale  nature  of  the 
physical  problem,  it  is  not  necessary  to  solve  all  four  equations  simultaneously  in  time.  Instead, 
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here  we  adopt  a  different  solution  strategy,  which  takes  place  in  an  iterative  manner  with  each 
iteration  containing  two  solution  phases. 

For  the  first  phase,  Maxwell’s  equations  together  with  the  momentum  transfer  equation  are 
solved  for  one  full  period  T  with  a  time  step  size  AtEM  to  obtain  the  root-mean-square  (RMS) 
value  of  the  electric  field  at  every  point  r  in  the  solution  domain 

/i  rto+T  2  \1/2 

Ems(r)  =  (^  J  \\E(r,t)  +  Emc(r,  t)||  dtj  (28) 

which  is  then  used  to  calculate  the  effective  electric  field  (with  o>  being  the  angular  frequency  of 
the  incident  field) 


V 1  +  W2/Um 


(29) 


and  the  reduced  electric  field  Eeff/p.  The  reduced  electric  field  is  finally  used  to  calculate  the 
plasma  parameters  such  as  the  ionization  and  attachment  frequency  u;  and  va,  and  the  plasma 
diffusion  coefficient  Deff  introduced  in  the  preceding  section,  which  are  all  functions  of  space 
and  time,  apparently. 


Once  all  the  plasma  parameters  are  obtained  at  every  point  in  the  solution  domain,  the  electron 
diffusion  equation  is  solved  for  one  time  step  with  a  time  step  size  AtPL  =  T,  which  is  the  second 
phase  in  the  iterative  scheme.  The  newly  obtained  electron  density  n  is  then  used  in  the 
momentum  transfer  equation  for  the  first  solution  phase  to  continue  for  another  period.  The 
entire  solution  thus  march  on  in  time  in  an  iterative  manner. 


It  is  worth  mentioning  that  to  make  the  entire  numerical  scheme  uniformly  high  order  accurate; 
the  time  integration  scheme  must  also  be  high  order  accurate.  To  this  end,  the  high-order  explicit 
Runge-Kutta  method  is  used  in  this  work. 

3.3  Mathematical  Basics  and  Considerations 


In  this  section,  the  mathematical  formulation  of  the  physical  problem  is  described,  with  related 
mathematic  basics  and  considerations  given.  To  achieve  the  flexibility  of  geometric  modeling, 
the  high-order  accuracy  in  both  spatial  and  temporal  discretization’s,  and  more  importantly,  the 
continuity  in  all  electromagnetic  field  components,  the  nodal  DGTD  method  is  employed  to 
solve  all  four  system  equations  given  in  the  preceding  section. 

3.3.1  Nodal  Discontinuous-Galerkin  Time-Domain  Method 

The  nodal  DGTD  method  used  in  this  research  work  is  different  from  the  vector  DG  method 
used  in  computational  electromagnetics  community  in  a  sense  that  the  physical  quantities  are 
defined  on  each  node  (will  be  discussed  in  a  following  section),  with  the  vector  quantities 
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expanded  into  x,  y,  and  z  components.  The  DGTD  used  in  this  work  for  Maxwell’s  equations  and 
for  electron  diffusion  equation  are  described  below. 

3.3. 1.1  Nodal  DGTD  for  the  Hyperbolic  PDE  -  Maxwell’s  Equations 

As  discussed  before,  Maxwell’s  equations  need  to  be  solved  along  with  the  momentum  transfer 
equation  for  an  entire  period  before  the  electron  diffusion  equation  can  be  solved.  Since  the 
momentum  transfer  equation  is  valid  point-wisely,  there  is  no  global  coupling  (information 
communication)  involved  in  the  solution.  Here  the  discussion  is  focused  on  the  solution  of 
Maxwell’s  equations. 

Consider  Maxwell’s  equations 


1 

H  = - V  x  E  (30) 

/t 

E  =  —(y  x  H  —  /c)  (31) 

subject  to  the  following  boundary  conditions 

nx  E  =  0,  n-  H  =  0,  re  rPEC  (32) 

nx  H  =  0,  n  -  E  =  0,  re  Ipj^Q  (33) 

nxE+ZnxnxH  =  0,  nxH  —  YnxnxH  =  0,  re  fABC  (34) 


where  rPEC,  rPMC,  and  rABC  denote  the  boundaries  for  the  perfect  electric  conductor  (PEC),  the 
perfect  magnetic  conductor  (PMC),  and  the  absorbing  boundary  condition  (ABC),  respectively. 

Testing  the  above  equations  with  Lagrange  polynomials  lt  =  wlt  (vv  =  x,  y,  or  z)  on  each 
tetrahedron  element  Ve,  and  applying  the  divergence  theorem  twice  yield  the  strong  form 


If  H  dV 
IfEdV 


1 

1 

e 


[  IfVxEdV  +  l  lrnx(E*-E)dS 

Ve  'dVe 

I  IfVxH  dV  +  l  li  ■  n  x  (H*  -  H)  dS  -  f  lrJc 

've  hve  Jve 


dV 


(35) 

(36) 


where  E*  and  H*  are  the  numerical  fluxes  defined  on  the  boundary  faces  of  each  tetrahedron 
element.  Two  choices  of  these  numerical  fluxes  are  commonly  applied.  One  is  the  central  flux, 
and  the  other  is  the  upwind  flux. 


To  facilitate  the  description,  we  define  the  average  {■}  and  jump  [[■]  functions  as  follows 


(a)  =  2  ( a+  +  a  ),  [a]  =  a+  —  a 


(37) 
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where  the  superscripts  -  and  +  indicate  the  inside  and  the  outside  of  a  surface,  respectively.  Note 
that  the  above  definitions  apply  to  both  cases  where  a  is  a  scalar  or  a  vector  variable. 

•  Central  Flux 

E*  =  {. E },  H*  =  {H}  (38) 

•  Upwind  Flux 

E*  =-^[{YE}  +  ^nxm)  (39) 

H*  =^}{{ZH]~\nXm)  (40) 

Using  either  of  the  above  numerical  flux  and  expanding  E  and  H  in  terms  of  /ith  order  Lagrange 
polynomial  basis  functions 

Np  Np 

E  =  ^lj  (xExj  +  yEyj  +  zEzj)  ,  H  =  'Yjlj  (xHxj +yHyj +  zHzj)  (41) 

j=l  7=1 

where  Np  =  Ilf=i  (P  +  Q/dl  stands  for  the  number  of  DoFs  in  a  J-dimensional  element,  the 
strong  form  of  the  Maxwell’s  equations  can  be  converted  into  the  matrix  form  representation 


MeHx  =  -  i  (S*EZ  -  SzEy  +  MefFxE ) 

(42) 

MeHy  =  ~  (: S*EX  -  SXEZ  +  MefFE ) 

(43) 

MeHz  =  -- {SexEy  -  SeyEx  +  Af/F/) 

\l 

(44) 

MeEx  =  i  {SyHz  -  SezHy  +  Af/F"  -  MeJcx) 

(45) 

MeEy  =  i  (SezHx  -  SXHZ  +  MefFyH  -  Me]cy ) 

(46) 

MeEz  =  i  (< SexHy  -  SeyHx  +  M/Fzh  -  Me7cz) 

(47) 

where  the  flux  terms  are  given  by 

Fe  =  i  n  x  [£■]  (48) 

Fh  =^nx[Hl  (49) 

for  the  central  flux,  and 
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1 


(50) 


Fe  = 

fh  = 


Y+  +  Y~ 
1 


Z+  +Z- 


(Y+  n  x  [£■]  +  n  x  n  x  [tf]) 
(Z+nx[W]-nxnx  [£]) 


for  the  upwind  flux,  while 


(51) 


[M%-  =  [  lj  dV 

JVe 

K]y  =  |  'i  ^  dl/ 

[M/1,,  =  f  'f iS 

1  JdVe 

are  mass,  stiffness,  and  facial  mass  matrices  defined  in  each  tetrahedron  element  Ve  and  its 
surface  dVe,  respectively. 

To  solve  the  above  equations,  the  momentum  transfer  equation  needs  to  be  solved 
simultaneously.  Since  it  is  valid  in  a  point-wise  manner,  no  numerical  flux  is  needed.  Using  the 
similar  discretization  method,  the  momentum  transfer  equation  can  be  converted  into  a  set  of 
matrix  equations  (w  =  x,  y,  or  z) 


(52) 

(53) 

(54) 


j  —  _y  1  _j_  6  U  ( £  +  £inc)  (55) 

Jew  uvciJew  '  _  y-'w  '  cw  J  v  J 

m 

When  Maxwell’s  equations  and  the  momentum  transfer  equation  are  solved  together,  the  current 
values  of  the  physical  quantities  E,  H,  and  Jc  are  substituted  into  the  RHS  of  the  above  matrix 
equations  to  obtain  the  time  derivative  of  themselves,  which  are  then  used  to  obtain  their  updated 
values  at  the  next  time  step  using  an  explicit  time  integration  scheme. 

3.3.1.2  Local  DGTD  for  the  Parabolic  PDE  -  Electron  Diffusion  Equation 

The  diffusion  equation  has  a  very  different  mathematical  property  from  Maxwell’s  equations 
because  of  the  different  orders  in  the  spatial  derivative.  For  the  diffusion  equation,  if  a  similar 
method  as  described  in  the  preceding  section  is  adopted,  a  convergent  but  inconsistent  numerical 
solution  will  be  produced  [34],  Here,  the  term  “convergent”  is  used  to  refer  to  the  fact  that  when 
fixing  the  polynomial  order  but  increasing  the  spatial  discretization  density  (/z-refinement),  the 
numerical  result  will  converge  to  some  certain  solution;  the  term  “inconsistent”  is  used  to  refer  to 
the  fact  that  the  convergent  solutions  for  different  polynomial  orders  (/^-refinement)  are  different. 
As  a  result,  a  different  DG  scheme  needs  to  be  used  to  obtain  a  numerical  solution  that  is  both 
convergent  and  consistent. 

In  this  work,  the  local  discontinuous  Galerkin  (LDG)  method  in  time  domain  [17]  is  adopted, 
with  the  basic  idea  and  formulations  described  below. 
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Consider  the  diffusion  equation 


ii  -  V  ■  (DVu)  =  /  (r,  t)  e  fl  x  [0,  oo)  (56) 

subject  to  the  following  boundary  and  initial  conditions 

u{r,t)  =  gD  r  G  rD  (57) 

n  ■  DVu(r,  t)  =  qN  r  G  rN  (58) 

an  +  n  ■  DVu(r,  t)  =  gR  r  G  TR  (59) 

u(r,  0)  =  it0  r  G  H  (60) 


where  fl  denotes  the  three  dimensional  solution  domain,  rD,  rN,  and  TR  denote  the  Dirichlet, 
Neumann,  and  Robin  boundaries,  respectively.  To  solve  it  with  DG  method,  we  define  an 
auxiliary  variable  q  =  Vu  and  rewrite  the  above  equations  as  two  coupled  equations  with  only 
first  order  derivatives  in  space 


o 

II 

> 

1 

( r ,  t)  Gflx  [0,  oo) 

(61) 

u  —  V  ■  (Dq)  =  / 

(r,  t)  G  D  x  [0,  oo) 

(62) 

u(r,  t)  =  gD 

r  G  rD 

(63) 

n  ■  Dq(r,  t)  =  qN 

r  g  rN 

(64) 

au  +  fi  •  Dq(r,  t)  =  gR 

rer r 

(65) 

u(r,  0)  =  u0 

r  ell. 

(66) 

Testing  the  above  equations  with  Lagrange  polynomials  lt  and  =  wlt  (vv  =  x,  y,  or  z)  on  each 
tetrahedron  element  Ve,  and  applying  the  divergence  theorem  twice  yield  the  strong  form 


[  lrqdV=  [  lfVud.V-<f  lrn(u-u*)dS  (67) 

ve  've  hve 

[  lt  u  dV  =  [  ltV-  ( Dq )  dV-<fi  lt  n  ■  [Dq  -  ( Dq )*]  dS  +  f  kf  dV.  (68) 

Jve  Jve  hve  Jve 

In  the  above  strong  form,  u*  and  (Dq)*  are  the  numerical  fluxes  defined  on  the  boundary  faces 
of  each  tetrahedron  element.  Two  choices  of  these  numerical  fluxes  are  commonly  applied.  They 
are  given  below  using  the  average  {■}  and  jump  [■]  functions  defined  in  the  preceding  section. 

•  Central  Flux 


u *  =  {u},  (Dq)*  =  {Dq} 


•  LDG  Flux 


u*  =  {u}  +  ye  ■  ii  luj 
(Dq)*  =  {Dq}  +  pe  n-  [Dq]  -rje  n  [u] 


(69) 


(70) 

(71) 


where 
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npe  =  ^sign(n-r), 


1 

ny=  --  sign(n  ■  r) 


(72) 


with  r  being  a  given  constant  vector  to  mimic  the  wind  direction,  and  r\e  >  0  to  serve  as  a 
penalty  factor. 

Expanding  u  and  q  in  terms  of pth  order  Lagrange  polynomial  basis  functions 

Np  Np 

u  ^  '  lj  Uj  ,  q  ^  '  lj  (xqxj  +  yqyj  +  zqZj )  ( "73 ) 

j=l  7=1 

and  substituting  them  into  the  strong  form  to  obtain  the  matrix  form  representation  of  the 
diffusion  equation  ( w  =  x,  y,  or  z) 

Meqw  =  S^u  -  Mf[nw(u  -  u*)]  (74) 

Meu  =  ^  S^Dqw  -  Mf  ^  nw[Dqw  -  (J DqwT ]  +  Mef.  (75) 

W  W 

To  solve  the  diffusion  equation,  the  auxiliary  variable  q  is  first  calculated  using  the  information 
of  u,  and  the  time  derivative  it  is  then  obtained  using  the  values  of  both  u  and  q.  The  value  of  u 
at  the  next  time  step  can  be  obtained  with  a  time  integration  scheme  which  will  be  discussed 
later. 

3.3.2  Higher-Order  Spatial  Discretization 

After  the  discussion  on  the  formulation  for  different  equations,  the  high  order  basis  functions 
themselves  need  to  be  introduced.  In  this  section,  the  commonly  used  high  order  basis  functions, 
the  high  order  hierarchical  basis  functions  and  the  high  order  interpolatory  basis  functions,  are 
introduced,  and  both  basis  functions  are  used  in  the  development  of  our  coupled  DGTD  method. 
Since  tetrahedron  elements  are  used  in  the  geometric  discretization  of  the  solution  domain,  basis 
functions  that  are  defined  on  tetrahedrons  will  be  introduced.  Nevertheless,  those  defined  on 
different  types  of  elements  can  also  be  defined  in  a  similar  fashion. 

3.3.2.1  Definition  of  Basis  Functions:  Modal  Basis  and  Nodal  Basis 

Unlike  that  of  the  high-order  vector  basis  functions  widely  used  by  the  computational 
electromagnetics  community,  the  definition  of  the  high-order  scalar  basis  functions  is  a  bit 
different.  The  high-order  hierarchical  basis  functions,  which  are  more  commonly  known  as  the 
modal  basis,  are  a  set  of  orthonormal  polynomials  defined  on  the  simplex,  which  is  a  tetrahedron 
in  3D.  Specifically,  the  following  orthonormal  polynomial  of  order  p  >  i  +  j  +  k  >  0  are 
defined  using  the  well-known  Jacobi  polynomials 

<P„„Cr,s,  t)  =  V8P10'0(a)P/'+1'°<7>)(l  -  f>yPlt2‘+2,+2'<V)(l  -  c)1*’  (76) 
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where  a  =  —2  —  1,  b  =  2  —  —  1,  and  c  =  t,  with  the  3D  simplex  defined  by  r,s,t  >  — 1, 

r  +  s  +  t  <  —  1.  The  Jacobi  polynomials  are  given  by 


Pn,/?0)  =  ^2*n\  C1  ~  Z)  “(1  +  z)  ^^r{(1-z)a(1  +  z)^(1-z2)”}- 

When  a  =  /?  =  0,  the  Jacobi  polynomial  reduces  to  the  Legendre  polynomials 


dn 


(77) 


(78) 


Note  that  the  modal  basis  are  hierarchical  in  nature,  meaning  that  a  set  of  higher  order  modal 
basis  contains  all  of  the  lower  order  modal  basis  (similar  to  the  Fourier  basis). 


The  physical  quantity  can  therefore  be  expanded  in  terms  of  the  modal  basis  as  (here  the 
subscript  i,j,  k  are  lumped  into  n  for  simplicity) 

Np 

u(x,t)  =  I  Un(t)  Ipnix).  (79) 

n=l 

The  definition  of  high-order  interpolatory  basis  functions  has  multiple  choices.  For  example,  one 
can  use  monomials 


t,j,  k>0,i+j  +  k<p  (80) 

or  polynomials  as  the  basis  functions.  Here  the  widely  used  Lagrange  polynomials  are  used  as 
the  high-order  interpolatory  basis  functions.  One  unique  property  of  the  Lagrange  polynomial  is 
that  its  value  is  one  only  at  one  interpolating  node  (which  is  known  as  its  definition  node),  but 
zeroes  on  all  other  Np  —  1  interpolating  nodes.  It  is  for  this  reason  that  the  Lagrange  basis 
functions  are  better  known  as  nodal  basis,  as  opposed  to  the  modal  basis  introduced  earlier.  This 
property  has  two  pleasing  consequences.  One  is  that  the  expansion  coefficient  of  each  Lagrange 
basis  function  stands  for  the  value  of  the  physical  quantity  at  the  definition  node  of  the  basis 
function.  The  other  is  that  when  using  the  Lobatto  nodes  to  define  the  Lagrange  basis  functions 
(will  be  introduced  in  the  next  section),  the  value  of  a  Lagrange  basis  function  defined  in  a  d- 
dimensional  simplex  is  purely  zero  on  any  of  its  lower  d'-dimensional  simplex  subset  (0  <  d'  < 
d),  if  the  definition  node  of  the  basis  function  is  not  in  that  d' -dimensional  simplex.  For 
example,  for  the  Lagrange  basis  defined  inside  the  volume  of  a  tetrahedron,  its  value  is  purely 
zero  on  all  the  four  faces,  six  edges,  and  four  vertices  of  the  tetrahedron.  This  property  is 
especially  useful  when  implementing  the  DG  method,  because  when  trying  to  collect  the  values 
for  the  numerical  flux,  all  matters  are  the  DoFs  that  are  defined  on  the  common  triangle  of  two 
neighboring  tetrahedrons,  and  all  other  DoFs  have  zero  contribution. 

The  difficulty  of  applying  the  Lagrange  basis  functions  in  a  high  order  polynomial  expansion  is 
that  there  are  no  general  explicit  expressions  for  high  order  Lagrange  polynomials.  In  a 
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calculation  that  involves  the  Lagrange  polynomial,  what  one  can  do  is  to  relate  the  Lagrange 
basis  function  (nodal  basis)  to  the  orthonormal  basis  function  (modal  basis)  as 


Np  Np 

u(xu  t)  li(x)  =  ^  u„(t)  lpn(x)  (81) 

i= 1  71—1 

where  denotes  the  3D  Lagrange  polynomials  defined  on  the  interpolating  node  xt.  Matching 
the  values  of  u  at  all  the  interpolating  nodes  xt  (i  =  1, ... ,  Np)  yields 


u(jx,  t)  =  y 


Vu  =  u 


(82) 


where  Vnm  =  ipm(xn)  is  the  generalized  Vandermonde  matrix,  u  =  lit!, ... ,  uWp  j  ,  and  u  = 
|ux, ... ,  uWp|  .  The  Lagrange  polynomials  can  thus  be  expressed  as 

Np 

ln(x )  =  ^  (T2_1)mn  ipm(x)  ,  or  l  =  V~T  ip  (83) 

771=1 

and  can  be  used  to  calculate  the  elemental  matrices  which  will  be  discussed  in  a  following 
section. 

3.3.2.2  Choice  of  Interpolating  Nodes 

One  very  important  aspect  of  applying  the  Lagrange  basis  functions  is  how  to  choose  the 
interpolating  nodes.  If  a  set  of  equi-spaced  nodes  is  used  to  define  the  Lagrange  basis  functions, 
the  famous  Runge  phenomenon  will  lead  to  a  highly  oscillating  function  value  when  very  high 
order  polynomials  are  used.  The  oscillation  becomes  even  worse  near  the  boundary  of  the 
definition  domain  of  the  polynomial,  which  would  completely  ruin  the  DG  method  where  the 
numerical  flux  plays  a  central  role  of  the  algorithm.  To  remove  the  Runge  phenomenon,  a  set  of 
nonuniformly  spaced  nodes  is  usually  used  to  define  the  Lagrange  polynomial,  which  is  often  the 
root  set  of  certain  polynomials  defined  in  the  same  simplex.  Examples  of  such  polynomials 
include  the  Legendre  polynomial,  the  Chebyshev  polynomials  of  the  first  or  the  second  kind,  and 
the  Lobatto  polynomial.  For  the  application  of  DG  method,  the  Lobatto  nodes  are  especially 
attractive,  since  the  node  set  includes  those  on  the  boundary  of  a  simplex,  which  becomes  handy 
when  calculating  the  numerical  flux.  Unfortunately,  the  Lobatto  nodes  in  3D  simplex  are  not 
well-defined.  There  are  several  approaches  of  getting  the  Lobatto  nodes.  For  example,  to 
minimize  the  Lebesgue  constant,  this  is  a  quantity  that  measures  the  interpolation  quality  of  a 
given  polynomial,  or  to  maximizing  the  magnitude  of  the  determinant  of  the  Vandermonde 
matrix  using  either  constraint  optimization  method  or  Monte  Carlo  method.  In  this  work,  a  set  of 
optimized  Lobatto  nodes  in  3D  simplex  is  used,  as  shown  in  Figure  2  adopted  from  [19],  The 
optimized  set  of  nodes  produces  a  small  Lebesgue  constant  when  compared  to  other  node  sets 
available  in  the  literature,  and  reduces  to  ID  Gauss-Lobatto-Legendre  (GLL)  node  set  on  the 
edges  of  the  tetrahedron. 
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(d)  (e)  (f) 

Figure  2.  Optimized  Interpolating  Nodes  in  Tetrahedron  [19].  (a)  p  =  1,  (b)  p  =  2,  (c) 

p  =  4,  (d)  p  =  6,  (e)  p  =  8,  and  (f)  p  =  10 


Multiple  advantages  of  choosing  an  optimal  set  of  interpolating  nodes  can  be  achieved.  The  first 
one  is  that  it  results  in  a  smaller  Lebesgue  constant,  and  therefore,  a  better  interpolation  quality. 
Second,  it  also  reduces  the  spectral  radius  of  the  facial  matrix  M_17y,  where  7y  is  the  same  as 
Mf,  with  its  number  of  row  expanded  to  Np  to  allow  the  matrix-matrix  product  with  M-1.  Since 
the  spectral  radius  p(M_1T^)  determines  the  time  step  size  required  to  maintain  stability,  a 
smaller  spectral  radius  results  in  a  larger  time  step  size,  and  consequently,  a  more  efficient 
algorithm.  Shown  in  Figure  3  is  the  comparison  made  between  the  widely-used  hierarchical 
vector  basis  functions  and  the  orthonormal  basis  functions  used  in  this  work.  From  Figure  3(a), 
the  spectral  radius  of  the  orthonormal  basis  functions  is  smaller  than  that  of  the  vector  basis,  and 
hence  results  a  larger  step  size  as  can  be  seen  in  Figure  3(b). 
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(a)  (b) 

Figure  3.  (a)  Spectral  Radius  of  the  Facial  Matrix  Versus  Polynomial  Basis  Order;  (b) 
Time  Step  Size  Required  to  Maintain  Stability  Versus  Polynomial  Basis  Order.  The  Time 
Step  Sizes  Are  Normalized  Against  That  of  the  First-Order  Orthonormal  Basis  Function 


3.3.2.3  Elemental  Matrices  in  the  DG  Method 

After  the  high-order  basis  functions  and  the  associated  definition  nodes  are  obtained,  the 
elemental  matrices  in  the  DG  method  can  be  explicitly  expressed.  In  this  section,  the  calculation 
of  the  mass,  stiffness,  and  facial  mass  matrices  are  discussed. 

We  start  from  the  calculation  of  the  mass  matrix.  Using  the  relation  l  =  V~J  xf),  the  mass  matrix 
can  be  expressed  as 


Me  =Je  J  llT  dV=  Je  V~T  J  i/m/>t  dV  V"1  (84) 

where  I  is  the  reference  simplex  in  3D,  and  Je  is  the  Jacobian  matrix  mapping  from  the  reference 
simplex  to  the  real  tetrahedron,  and  is  reduced  to  a  single  scalar  if  the  geometric  discretization  is 
linear  (1st  order  tetrahedron).  Since  the  modal  basis  are  orthonormal,  the  mass  matrix  can  be 
written  in  an  explicit  form  as 


Me  =  Je  (WT)-1  =  Je  M.  (85) 

Note  that  in  the  explicit  form  of  the  mass  matrix,  no  numerical  quadrature  is  required.  The  facial 
mass  matrix  can  be  obtained  in  a  similar  manner;  the  only  difference  is  that  everything  in  the 
explicit  expression  reduces  from  3D  to  2D  as 

Mf  =  ;;  (VfVjy1  =  ;;  uf.  (86) 
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The  stiffness  matrix  involves  the  partial  derivative  with  respect  to  x,  y,  or  z.  Here  we  take  Sx  for 
example  to  illustrate  its  calculation. 


dlj  f  fdljdr  dhds  dh  dt\ 

li  yr~  dV  =  I  /;  (  -t—  — — I-  -t—  — — I-  -r— —  I 
dx  \dr  dx  ds  dx  dtdxj 


dV 


=  rv 


l 


dlj  r  dlj  f  dlj 


(87) 


In  the  above  derivation,  the  derivatives  rx,  sx,  and  tx  are  constants  if  linear  tetrahedron  elements 
are  used,  and  the  integration  terms  can  be  further  written  as 


1 


dl, 

h  In  dV 


^ '  Mjn  Drjlj 

n= 1 


(88) 


where  the  interpolation  of  using  Lagrange  polynomials  is  applied,  and  Dr  nj  is  the  differential 
matrix  with  respect  to  the  variable  r 


D 


r,nj 


(89) 


which  can  also  be  expressed  analytically  using  the  relation  between  the  Lagrange  and  the 
orthonormal  basis  functions.  Details  will  not  be  given  here.  The  stiffness  matrix  can  finally  be 
expressed  as 


Sx  —  M  (rxDr  +  sxDs  +  txDt )  (90) 

and  Sy  and  5|  can  be  expressed  similarly. 

3.3.3  Higher  Order  Temporal  Discretization 

After  the  system  equations  are  discretized  with  high  order  basis  functions  in  space,  a  set  of  ODEs 
can  be  obtained  and  formally  written  as 

ii  =  F(t,u ),  u(t0)  =  uQ.  (91) 

The  system  of  ODEs  can  be  solved  using  the  well-known  explicit  high  order  Runge-Kutta 
methods  [35],  for  instance,  the  classic  fourth  order  Runge-Kutta  method,  to  achieve  the  high 
order  accuracy  in  time.  However,  the  traditional  Runge-Kutta  methods  are  in  general  not  strong- 
stability  preserving  (SSP),  and  the  SSP  coefficient  (which  is  the  ratio  between  the  maximum 
time  step  size  for  a  given  time  integrator  and  that  for  the  forward  Euler’s  method)  C  *  1. 

In  this  work  we  employ  the  SSP  Runge-Kutta  (SSPRK)  method  [36],  which  is  also  known  as  the 
total  variation  diminishing  (TVD)  Runge-Kutta  Method  [3  7]  [3  8]  in  the  past.  The  strong-stability 
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preserving  property  means  that  if  the  embedded  first  order  Euler’s  method  is  stable,  the  high 
order  SSPRK  method  is  also  stable,  which  diminishes  the  total  variation  in  the  solution.  In  the 
meantime,  the  SSP  coefficient  can  be  made  larger  to  allow  a  larger  time  step  size.  For  example, 
the  low-storage  four-stage  third-order  SSPRK  is  given  by 


<7i 

<72 

<?2 

<72 


=  un 


h 


=  <7i  +  2F(<7i) 

=  <72  +  2F^2) 


-3^1  +  3 

=  <72  +^F(<72) 


h 


<72  +  2F(<?2) 


(92) 

(93) 

(94) 

(95) 

(96) 


and  the  SSP  coefficient  C  =  2.  In  this  work,  a 
employed  because  of  its  large  SSP  coefficient  C 
SSPRK  methods. 


low-storage  ten-stage  fourth-order  SSPRK  is 
=  6  which  makes  it  more  efficient  than  other 


Despite  its  attractive  properties,  the  SSPRK  methods  have  their  own  issue  that  they  are  originally 
designed  for  PDE  with  a  constant  coefficient  operator  F(u),  which  means  on  the  RHS  of  the 
equation,  no  time  variable  can  be  involved.  When  looking  back  to  the  coupled  system  equations 
we  are  dealing  with,  it  is  realized  that  the  SSPRK  methods  can  be  applied  directly  to  Maxwell’s 
equations  and  the  electron  diffusion  equation.  However,  since  the  incident  field  directly  involves 
the  time  variable  in  its  definition,  SSPRK  cannot  be  applied  in  a  straightforward  way.  In  this 
work,  this  issue  is  resolved  by  defining  auxiliary  variables  as  detailed  next. 

Since  Emc  =  E 0  sin[oi(t  —  r)]  =  i;0[cos(a)T)  sin(o)t)  —  sin(o)r)  cos(o)t)],  where  r  is  the  time 
delay  associated  with  each  spatial  point  r  with  respect  to  a  reference  point  r0, 
r  =  (r  —  r0)  ■  kmc/c,  by  defining  two  auxiliary  variables  =  sin(oit)  and  f2  =  cos  (cat),  a 


new  system  with  a  constant  coefficient  operator  is  obtained  as 
1 

H  = - V  x  E  (97) 

F 

E  =  |(Vxff-/c)  (98) 

2  2 

/c  =  -vmJc  +  —E  +  —E0[cos(mt)  A  -  sin(o)r)  f2\  (99) 

m  m 

A  =  cof2  (100) 

A  =  -toft.  (ioi) 


Therefore,  the  SSPRK  method  can  be  applied  and  the  high-order  accuracy  in  both  time  and  space 
is  achieved  with  a  relaxed  stability  condition. 
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3.4  Computational  Acceleration  Techniques 


As  pointed  out  in  the  preceding  section,  the  physical  problem  we  are  dealing  here  is  a  multiscale 
problem,  which  requires  an  extremely  intensive  computation  to  complete  the  whole  simulation. 
In  order  to  make  the  method  practical,  two  very  effective  acceleration  techniques  have  been 
developed  in  this  work.  One  uses  the  software-based  (algorithm-based)  acceleration  technique; 
the  other  uses  the  hardware-based  acceleration  technique. 

3.4.1  Dynamic  /^-Adaption 

To  resolve  the  fast  varying  physical  quantity,  one  approach  is  to  use  higher-order  basis  functions, 
which  is  known  as  the  /^-refinement.  However,  the  application  of  a  uniformly  high-order  basis 
functions  will  result  in  a  very  large  number  of  DoFs  in  the  system,  which  makes  the  simulation 
prohibitively  long.  In  the  simulation  of  the  time  domain  evolution  of  the  plasma,  the  plasma 
front,  where  the  largest  density  gradient  is  located,  needs  the  highest  spatial  resolution,  and  is 
moving  with  time,  while  the  rest  of  the  plasma  region  is  relatively  smooth  and  slowly  varying. 
To  take  the  advantage  of  the  higher-order  basis  function  and  control  the  overall  computational 
cost,  a  flexible  dynamic  /> adaption  algorithm  is  developed  in  this  work.  Instead  of  using  a 
uniform  polynomial  order  throughout  the  entire  solution  domain,  the  algorithm  first  identifies 
regions  that  need  higher  spatial  resolution,  (and  therefore,  higher  polynomial  order)  versus  those 
do  not  need  such  high  order  spatial  resolution.  Since  the  plasma  is  evolving  and  propagating  in 
space  and  time,  this  identification  process  is  made  dynamic  such  that  the  order  of  a  given 
tetrahedron  element  is  adjusted  according  to  the  specific  need. 

The  important  point  in  the  dynamic  adjustment  is  to  find  a  good  measuring  quantity  to  determine 
the  polynomial  order.  Take  the  electron  density  for  example,  since  the  high  order  polynomial  is 
employed  to  resolve  its  spatial  variation,  how  fast  it  actually  varies  in  space  would  be  a  good 
indicator  of  the  polynomial  order  needed.  As  a  result,  the  magnitude  of  the  gradient  of  the 
electron  density  turns  out  to  be  the  intended  indicator  for  the  dynamic  /^-adaption  algorithm. 
Shown  in  Figure  4  is  the  contour  plot  of  the  electron  plasma  density  and  its  gradient,  which 
reveals  the  fact  that  around  the  plasma  front,  the  electron  density  decreases  fastest  and  therefore 
need  a  higher  order  polynomial  for  it  to  be  resolved. 

In  addition  to  determining  the  polynomial  order  based  on  the  density  gradient,  other  constraints 
can  also  be  set.  For  example,  one  can  require  that  the  order  difference  between  two  adjacent 
tetrahedrons  be  smaller  than  some  preset  value  to  avoid  sudden  change  of  spatial  resolution  and 
the  resulting  numerical  instability. 

The  process  of  dynamic  p- adaption  is  summarized  as  following  steps. 

1.  Set  initial  conditions; 

2.  Dynamic  /^-adaption:  change  basis  orders  adaptively  based  on  ||Vn||  and  other 
constraints; 

3.  Elevate  elemental  DoFs:  n,  E,  H,  and  Jc ; 

4.  Change  facial  connectivity  array; 
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5.  Solve  the  system  equations:  elevate  facial  DoFs  at  every  stage  of  the  SSPRK  method  to 
calculate  numerical  fluxes; 

6.  Go  to  step  2,  march  on  in  time  until  the  problem  is  solved. 


Figure  4.  Contour  Plot  of  the  Electron  Density  (Left  Half)  and  it’s  Gradient  (Right  Half) 


gradN 


le+24 


3.4.2  GPU  Acceleration 

The  DGTD  method  is  in  fact  an  element  level  domain  decomposition  method,  which 
decomposes  the  entire  solution  domain  into  discretization  elements,  and  solves  the  governing 
equations  within  each  of  the  elements.  With  the  aid  of  the  explicit  Runge-Kutta  method,  all  the 
information  needed  from  the  neighboring  elements  is  readily  obtained  from  the  previous  time 
steps.  As  a  result,  the  DGTD  method  can  be  parallelized  with  a  very  high  efficiency.  In  this 
work,  the  coupled  DGTD  algorithm  is  parallelized  using  graphics  processing  units  (GPUs)  using 
NVIDIA’s  CUDA  (compute  united  device  architecture)  parallel  programming  models  [39]. 

3.4.2.1  Single  GPU  Acceleration 

A  big  advantage  for  the  DGTD  method  to  be  accelerated  by  GPU  is  that  it  does  not  require  the 
storage  of  elemental  matrices.  As  described  in  Section  IV-2.3,  the  matrices  that  need  to  be  stored 
are  those  defined  on  the  reference  simplex  M,  Mf,  Dr,  Ds,  and  Dt.  The  element-related  variables 
are  only  several  constants  for  each  element,  which  include  Je,  Jj,  rxyz,  sxyz,  and  txy  z.  This 
nice  property  of  the  DGTD  method  (for  linear  elements)  makes  it  extremely  storage-efficient.  As 
a  result,  even  a  very  large  problem  can  be  solved  on  a  single  GPU  card  where  the  memory  is 
limited.  The  implementation  of  the  nodal  DGTD  method  on  GPU  is  similar  to  those  reported  in 
[40] [41],  and  will  not  be  detailed  here.  When  the  GPU  is  enabled,  the  computational  speed-up  is 
about  10  times  compared  to  the  OpenMP -parallelized  CPU  code  running  on  eight  CPU  cores. 
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3.4.2.2  GPU  Acceleration  of  the  DGTD  Algorithm  with  /^-Adaption 

Since  the  dynamic  /(-adaption  algorithm  is  developed  in  CPU-based  code,  it  is  desired  to  extend 
this  capability  to  GPUs  for  a  better  efficiency.  However,  the  extension  from  a  GPU  code  with 
uniform  polynomial  order  to  dynamic  orders  is  quite  challenging,  because  the  flexibility  of  the  p- 
adapted  code  would  prevent  the  GPU  implementation  from  having  a  high  parallel  efficiency.  In 
this  work,  the  GPU  algorithm  with  dynamic  p- adaption  technique  is  carefully  designed  to 
maximize  its  parallel  efficiency.  The  matrices  for  all  possible  polynomial  orders  that  are  defined 
on  the  reference  simplex  are  pre-stored  in  the  GPU  global  memory  to  allow  the  potential  matrix- 
vector  products  for  the  corresponding  order.  In  addition,  the  elevation  matrices  that  elevate 
elemental  DoFs  between  different  orders  are  also  pre-stored  in  the  GPU  global  memory.  At  a 
specific  time  step  in  the  time-stepping  process,  the  entire  solution  domain  is  first  scanned,  and 
the  DoFs  in  different  tetrahedrons  are  elevated  (either  increased  or  decreased)  to  the  correct 
polynomial  order  whenever  and  wherever  needed.  The  numerical  fluxes  on  the  surfaces  that 
interface  with  tetrahedrons  of  different  polynomial  orders  are  then  calculated  to  permit  the 
correct  information  exchange  between  neighboring  tetrahedrons.  The  tetrahedrons  with  the  same 
polynomial  order  are  then  aligned  in  the  same  streaming  array  and  processed  together  since  they 
are  all  multiplied  by  the  same  matrices,  to  obtain  an  optimized  speed-up.  The  dynamic  /(-adapted 
GPU  algorithm  can  eventually  achieve  about  eight-time  speed-up  compared  with  the  dynamic  p- 
adapted  CPU  algorithm  running  in  parallel  on  eight  CPU  cores. 

4.0  RESULTS  AND  DISCUSSION 

Numerical  results  are  given  in  this  section  to  demonstrate  the  accuracy  of  the  DGTD  method  and 
model  several  physical  phenomena  related  to  the  high-power  microwave  breakdown  and  the 
plasma  formation  and  evolution.  The  extremely  high  gradient  of  the  electron  density  distribution 
is  resolved  by  using  either  a  dense  geometrical  discretization  or  a  set  of  higher  order  basis 
functions. 

4.1  Validation  of  the  LDG  Method  for  the  Diffusion  Equation 

The  first  example  given  here  is  the  solution  of  the  diffusion  equation  using  the  LDG  method  for 
the  validation  purpose.  Consider  the  following  density  diffusion  problem  with  a  diffusion 
coefficient  D  =  1 


ii  —  V  ■  (DVu)  =  0 

(r,  t)  E  fix  [0,  oo) 

(102) 

u(r,  t)  =  0 

r  G  rD 

(103) 

u(r,  0)  =  u0 

r  E  fl 

(104) 

defined  in  a  spherical  solution  domain  fl  where  r  G  [0,  a]  (a  =  2  mm).  The  initial  condition  of 
the  density  distribution  is  a  Gaussian  function  with  a  standard  deviation  a  =  0.4  mm 


u0  =  e  2cr2 . 


(105) 


This  problem  has  an  analytical  solution  given  by 
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(106) 


u(r,t)  =  ^  cnj0(pnr )  e 

n= 1 

where  j0  is  the  zeroth  order  spherical  Bessel  function  of  the  first  kind,  /?n  =  — ,  and  cn  can  be 
found  by  matching  the  initial  condition. 

This  problem  is  solved  numerically  using  the  LDG  method  in  time  for  a  total  of  Nt  =  13150 
time  steps  with  a  time  step  size  At  =  20  ps.  The  final  result  along  the  center  line  of  the  solution 
domain  at  t  =  263  ns  is  compared  with  the  analytical  solution  and  is  shown  in  Figure  5,  from 
which  excellent  agreement  can  be  observed.  Figure  6  shows  the  contour  plot  of  the  equal-density 
surfaces  of  the  electron  density  distribution  at  the  final  time  step.  It  can  be  seen  from  this  figure 
that  even  after  time  stepping  for  over  13  thousand  steps,  the  equal-density  surfaces  of  the  density 
distribution  given  by  the  numerical  solution  is  still  very  smooth,  which  demonstrates  the  high 
accuracy  of  the  LDG  method  employed  in  this  project. 


Coordinate  (m) 

Figure  5.  Comparison  between  the  Numerical  and  the  Analytical  Solutions  of  the  Diffusion 

Equation 
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Figure  6.  Contour  Plot  of  the  Density  Distribution  at  the  Final  Time  Step 


4.2  Modeling  of  Plasma  Shielding  of  the  Electromagnetic  Wave 

In  this  example,  the  air  breakdown  under  the  high-power  microwave  operation  and  the  formation 
and  evolution  of  the  plasma  is  simulated.  This  problem  is  adopted  from  [33]  and  was  originally 
simulated  with  FDTD  and  FVTD  (finite  volume  time  domain)  methods  in  two  dimensions.  In 
this  work,  this  example  is  reproduced  using  the  DGTD  method  in  three  dimensions. 

4.2.1  Physics  Involved 

During  the  high-power  microwave  operation,  the  electric  field  with  a  high  power  density  will 
ionize  neutral  gas  and  produce  electrons  and  positive  ions.  If  the  power  density  of  the  electric 
field  is  high  enough,  such  an  ionization  process  will  continue,  which  results  in  an  exponential 
increase  of  the  electron  density  known  as  breakdown.  However,  this  process  will  not  go  on 
forever.  According  to  the  basic  plasma  physics,  when  the  plasma  density  increases  to  a  point  that 
the  corresponding  plasma  frequency  becomes  greater  than  the  incident  frequency,  the  electrons 
in  the  plasma  bulk  can  oscillate  fast  enough  such  that  the  secondary  fields  radiated  from  the 
electron  oscillation  cancels  the  primary  (incident)  fields.  Macroscopically,  the  plasma  bulk  is 
acting  like  a  piece  of  electromagnetic  shield  which  prevents  the  electromagnetic  wave  from 
propagating  through.  This  leads  to  a  decrease  of  the  total  electric  field  intensity  in  the  plasma 
bulk  which  lowers  the  ionization  rate.  As  a  result,  the  electron  density  will  get  saturated  at  a 
certain  level  and  the  entire  system  will  run  into  a  steady  state. 

4.2.2  Numerical  Settings 

As  shown  in  Figure  7,  the  solution  domain  considered  in  this  example  is  a  parallel  plate 
waveguide  with  a  metallic  wall  placed  in  between  of  the  two  parallel  plates,  which  forms  a 
rectangular  aperture  (denoted  in  red  in  the  figure).  The  solution  domain  is  truncated  from  the  left 
and  the  right  using  the  ABC.  Different  from  the  2D  problem  considered  in  [33],  in  this  work  a 
3D  problem  is  solved  by  setting  a  thickness  to  the  structure  (the  dimension  not  shown  in  Figure 
7)  and  placing  the  PMC  boundary  conditions  on  both  sides  in  the  thickness  direction  to 
reproduce  the  2D  phenomenon  in  3D.  A  set  of  25-GHz,  2.0-MV/m,  vertically  polarized  plane 
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wave  with  a  tapered  sinusoidal  temporal  profile  is  launched  from  the  left  boundary  and 
propagates  toward  the  right  direction  (the  z  direction). 


1.85  0.30  3.85 

I*  I*  *1 


Figure  7.  Illustration  of  the  Solution  Domain  of  the  Microwave  Breakdown  Problem  in  a 

Metallic  Aperture 


Two  cases  are  considered  and  compared  in  this  example.  One  is  the  vacuum  case  where  the 
entire  structure  is  assumed  to  be  filled  with  vacuum,  and  therefore,  the  metallic  wall  is  simply 
reflecting  a  portion  of  the  incident  wave  and  enhancing  the  fields  in  the  aperture  area.  The  other 
is  the  air  case  where  the  100-torr  low  pressure  air  is  assumed  to  be  confined  in  the  aperture  area, 
with  the  rest  of  the  domain  filled  with  vacuum.  When  the  plane  wave  is  incident  upon  the  air 
region,  it  is  enhanced  by  the  metallic  wall  and  triggers  the  air  breakdown.  In  the  numerical 
simulation,  four  boundaries  of  the  air  aperture  are  set  as  homogeneous  Dirichlet  boundaries  for 
the  electron  density,  and  the  other  two  boundaries  in  the  thickness  (x)  direction  are  set  as 
homogeneous  Neumann  boundaries.  To  satisfy  the  Dirichlet  boundary  condition,  the  initial 
spatial  profile  of  the  electron  density  is  set  to  be  sinusoidal  in  both  the  vertical  (y)  and  the 
horizontal  (z)  directions,  which  reaches  its  maximum  value  of  1015/m3  at  the  center  of  the 
aperture.  In  the  x  direction,  the  electron  density  distribution  is  uniform. 

In  both  cases,  the  electric  field  is  recorded  at  two  observation  points  PI  (x  =  y  =  0.0  mm, 
z  =  0.0  mm)  and  P2  (x  =  y  =  0.0  mm,  z  =  2.81  mm)  shown  in  Figure  7.  In  the  air  case,  the 
electron  density  n  and  the  effective  field  Etn  in  the  air  region  are  also  recorded  to  demonstrate 
the  spatial  and  temporal  evolution  of  the  plasma. 

4.2.3  Numerical  Results  and  Analysis 

Shown  in  Figure  8  is  the  temporal  profile  of  the  electric  fields  recorded  at  PI  and  P2  in  both  the 
vacuum  and  the  air  cases.  Different  behaviors  can  be  observed  in  the  two  cases.  In  the  vacuum 
case,  the  electric  field  intensity  at  both  observation  points  follows  the  tapered  sinusoidal  profile 
of  the  incident  wave,  with  the  field  at  PI  point  enhanced  by  the  narrow  aperture,  and  that  at  P2 
point  weakened  comparing  with  the  magnitude  of  the  incident  wave  (2  MV/m)  due  to  the  partial 
reflection  of  the  incident  wave  by  the  metallic  wall.  In  the  air  case,  however,  the  electric  field 
intensity  recorded  at  both  observation  points  increases  in  the  early  time,  but  decreases  suddenly 
at  PI  point  and  gradually  at  P2  point,  and  becomes  stabilized  at  much  lower  values  after  0.5  ns. 
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Time  (n&)  Time  (ns) 

(a)  (b) 

Figure  8.  Temporal  Response  of  the  Electric  Fields  Recorded  at  PI  and  P2  in  (a)  the 

Vacuum  and  (b)  the  Air  Cases 


The  electric  fields  recorded  in  both  cases  at  P2  point  are  compared  in  Figure  9(a),  from  which  it 
is  clear  that  the  electric  field  observed  in  the  air  case  agrees  well  with  that  in  the  vacuum  case  up 
to  0.35  ns.  After  that,  the  electric  field  in  the  air  case  starts  to  get  shielded  and  its  magnitude 
decreases  and  gets  stabilized  after  0.5  ns.  From  the  electron  density  recorded  at  PI  point  shown 
in  Figure  9(b),  it  can  be  seen  that  initially  the  electron  density  increases  exponentially  from  its 
initial  value  of  1015/m3  to  over  1021/m3  within  the  first  0.5  ns,  and  gets  saturated  afterwards. 


Time  (ns)  Time  (ns) 

(a)  (b) 

Figure  9.  (a)  Comparison  of  the  Temporal  Response  of  the  Electric  Fields  Recorded  at  P2 
in  Both  Cases;  (b)  Temporal  Evolution  of  the  Electron  Density  Recorded  at  PI  in  the  Air 

Case 


To  obtain  a  better  understanding  of  the  physical  process,  the  spatial  and  temporal  evolution  of 
the  electron  density  distribution  is  presented.  Shown  in  Figure  10  is  the  electron  density 
distribution  in  the  air  aperture  at  0.28  ns.  From  its  variation  in  the  vertical  direction  shown  in 
Figure  10(a),  two  peaks  near  the  boundaries  are  observed,  which  is  due  to  the  edge  singularity 
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that  enhances  the  electric  field  near  these  two  boundaries.  To  satisfy  the  Dirichlet  boundary 
condition  at  the  two  boundaries,  the  electron  density  has  to  have  very  large  density  gradients 
around  the  boundaries,  which  need  to  be  resolved  with  high  spatial  resolution.  From  Figure 
10(b),  the  electron  density  distribution  in  the  horizontal  direction  is  almost  sinusoidal  except  for 
a  small  shift  towards  the  right,  which  is  due  to  the  nonuniform  field  distribution  that  results  in 
different  ionization  rate  along  this  direction.  A  two  dimensional  snapshot  of  the  electron  density 
at  0.28  ns  is  given  in  Figure  10(c)  to  provide  a  whole  picture  of  the  density  distribution. 


(a)  (b)  (c) 

Figure  10.  Electron  Density  Distribution  in  the  Air  Aperture  at  0.28  ns.  (a)  One 
Dimensional  Distribution  Along  y  Direction;  (b)  One  Dimensional  Distribution  Along  z 

Direction;  (c)  Two  Dimensional  Snapshot 


Spatial  and  temporal  evolution  of  the  electron  density  n  and  the  effective  field  Ee ff  from  0.28  ns 
to  4.00  ns  is  presented  in  Figure  11  in  terms  of  two  dimensional  snapshots  once  every 
electromagnetic  period  (1/25  GHz  =  0.04  ns).  At  0.28  ns  [Figure  11(a)],  the  electron  density  is 
too  low  to  disturb  the  incident  field,  and  the  effective  field  reaches  its  maximum  at  four  comers 
(hot  spots)  of  the  air  aperture  due  to  the  edge  singularity.  These  four  hot  spots  lead  to  faster 
ionization  and  result  in  the  four  areas  with  higher  electron  density  that  can  be  observed  from 
Figure  11(b)  at  0.32  ns.  These  four  plasma  areas  in  return  begin  to  shield  the  electromagnetic 
field  within  the  areas  and  in  the  same  time  radiate  at  their  boundaries  secondary  fields  that 
enhance  the  incident  fields,  resulting  smaller  effective  field  within  the  plasma  bulk  and  new  hot 
spots  around  it.  The  new  hot  spots  ionize  more  air  that  forms  new  areas  of  plasma  bulk  which 
shields  more  electromagnetic  fields  and  generates  newer  hot  spots  [Figure  1 1(c)],  Such  a  process 
keeps  on  going  until  about  0.5  ns,  when  the  ionized  plasma  pattern  is  fully  developed  and  the 
electromagnetic  field  is  shielded  in  the  entire  air  aperture,  resulting  in  a  significant  magnitude 
drop  of  the  observed  electric  fields  at  both  PI  and  P2  points,  as  those  presented  in  Figure  8(b) 
and  Figure  9(a). 

After  0.5  ns,  the  ionization  of  the  neutral  air  becomes  very  slow.  Instead,  the  dominant  physical 
mechanism  changes  from  ionization  to  diffusion,  where  the  electron  diffuses  from  the  higher- 
density  regions  to  lower-density  regions.  Specifically,  the  plasma  front  exhibits  the  free  diffusion 
with  a  high  electron  diffusion  coefficient,  and  the  plasma  bulk  exhibits  the  ambipolar  diffusion 
which  has  a  much  lower  ambipolar  diffusion  coefficient.  As  a  result,  from  0.60  ns  [Figure  1  l(i)] 
to  4.00  ns  [Figure  ll(j)],  the  overall  pattern  of  the  plasma  formation  does  not  change 
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significantly  due  to  the  slow  ambipolar  diffusion.  But  the  boundaries  of  the  plasma  regions 
become  blurred  and  smoother,  because  of  the  free  electron  diffusion  near  the  plasma  boundaries 
that  inject  the  lower-density  regions  with  more  electrons. 
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(i)  0) 

Figure  11.  Snapshots  of  the  Time  Evolution  of  the  Electron  Density  and  the  Effective 
Electric  Field  in  the  Air  Aperture  at  (a)  0.28  ns,  (b)  0.32  ns,  (c)  0.36  ns,  (d)  0.40  ns,  (e)  0.44 
ns,  (f)  0.48  ns,  (g)  0.52  ns,  (h)  0.56  ns,  (i)  0.60  ns,  and  (j)  4.00  ns,  Respectively 
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From  this  example,  it  is  clear  that  the  methodology  developed  in  this  work  is  fully  capable  to 
capture  the  dominant  physical  mechanism  (ionization  versus  diffusion)  and  provide  deeper 
insight  to  the  physical  phenomenon  and  process. 

4.3  Reproduction  of  the  Microwave  Streamer 

In  the  last  example,  the  generation  of  the  microwave  streamer  is  considered.  As  shown  in  Figure 
12,  the  microwave  streamer  discharge  demonstrates  the  elongation  effect  of  a  plasma  bulk  under 
the  excitation  of  a  standing  electric  field.  Under  the  combined  effect  of  diffusion,  ionization, 
attachment,  and  recombination,  the  plasma  bulk  will  grow  longer  from  its  initial  shape,  for 
example,  a  Gaussian  distribution,  in  the  direction  of  the  externally  applied  electric  field.  The 
generation  of  the  microwave  streamer  is  experimentally  observed  and  reported  in  [42], 


Figure  12.  The  Microwave  Streamer  Discharge  Recorded  in  Experiments.  Streamer 
Discharges  in  Air  at  (a)  480  torr  and  (b)  760  torr.  Streamer  Discharges  in  Hydrogen  at  (c) 

480  torr  and  (d)  1000  torr  [42] 


4.3.1  Physics  Involved 

Shown  in  Figure  13  is  the  effective  ionization  frequency  and  electron  free  diffusion  coefficient  as 
functions  of  reduced  effective  field.  When  the  effective  field  is  below  the  critical  value  (shown 
as  the  vertical  dashed  line),  the  attachment  effect  dominates  over  the  ionization  effect,  while 
when  the  effective  field  is  higher  than  the  critical  value,  the  ionization  frequency  increases 
dramatically  and  becomes  dominant.  The  diffusion  coefficient  is  also  a  function  of  the  reduced 
effective  field,  but  it  increases  much  slower  compared  to  the  ionization  frequency. 
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Figure  13.  Effective  Ionization  Frequency  and  Electron  Diffusion  Coefficient  as  Functions 

of  Effective  Field 


4.3.2  Numerical  Settings 

In  the  experiment,  for  the  air  breakdown  to  begin,  a  highly  intensive  Gaussian  beam  is  shined 
into  an  air  chamber,  which  initiates  the  air  breakdown  at  its  focal  point.  In  the  numerical 
simulation,  the  simulation  domain  is  set  as  a  3D  cylindrical  region  filled  with  760-torr  air,  with 
all  its  boundaries  set  as  ABC  for  the  electromagnetic  waves  and  homogeneous  Dirichlet 
boundary  condition  for  the  plasma.  Instead  of  the  Gaussian  beam,  a  110-GHz,  6-MV/m, 
horizontally  polarized  standing  wave  is  launched  in  the  solution  domain  with  the  maximum 
value  of  the  standing  wave  located  along  the  center  line  (the  symmetric  axis)  of  the  cylindrical 
air  box,  and  a  seed  electron  bulk  is  placed  at  the  center  point  of  the  air  box,  which  has  a  density 
distribution  as  a  Gaussian  dot  with  the  maximum  density  set  as  1015/m3  and  the  standard 
deviation  set  as  100  jum,  as  shown  in  Figure  14. 
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Figure  14.  The  Initial  Condition  of  the  Electron  Distribution,  a  Gaussian  Dot 
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4.3.3  Numerical  Results  and  Analysis 


The  spatial  and  temporal  evolution  of  the  microwave  streamer  is  shown  in  Figure  15  at  multiple 
time  steps.  From  these  figures,  it  can  be  seen  that  from  the  initial  time  step  [Figure  15(a)]  to 
about  12.0  ns  [Figure  15(b)],  the  shape  of  the  plasma  bulk  remains  almost  the  same,  while  the 
electron  density  increases  from  the  initial  1015/m3  to  about  1021/m3  due  to  ionization  of  neutral 
gas  into  the  plasma.  Since  the  secondary  fields  radiated  by  the  plasma  oscillation  are  still  too 
weak  to  disturb  the  incident  waves,  the  effective  fields  shown  in  these  two  figures  are  almost 
straight  lines  with  similar  magnitude,  indicating  that  the  effective  field  is  dominated  by  the 
incident  standing  wave.  From  12.0  ns,  the  electron  density  grows  to  a  significant  level  such  that 
its  radiated  fields  start  to  interact  with  the  incident  waves,  which  results  in  the  increase  of  the 
effective  fields  at  the  plasma  front  in  the  direction  of  the  electric  field  (the  horizontal  direction  in 
the  figures)  and  the  decrease  in  between  of  these  two  peaks  of  the  effective  fields  due  to  the 
plasma  shielding  to  the  electromagnetic  fields,  as  can  be  observed  from  Figure  15(c)  at  14.0  ns. 
The  difference  of  the  effective  fields  in  the  plasma  region  result  in  different  ionization 
frequencies  and  diffusion  coefficients,  which  leads  to  a  faster  ionization  rate  and  diffusion 
velocity  at  the  two  tips  of  the  plasma  region  compared  with  the  inside  of  the  plasma  bulk, 
causing  the  electron  density  at  these  regions  to  increase  faster  than  that  at  other  regions.  This 
makes  the  plasma  bulk  change  its  shape  and  grow  longer  in  the  horizontal  direction.  Such  a 
process  continues  and  generates  even  sharper  peaks  of  the  effective  fields  at  the  plasma  front  and 
even  lower  magnitude  of  the  effective  fields  inside  the  plasma  bulk,  and  therefore,  the  elongation 
of  the  plasma  bulk  continues,  which  eventually  produces  the  so-called  microwave  streamer  in  the 
direction  of  the  incident  electric  field  as  can  be  seen  in  Figure  15(d). 
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Figure  15.  The  Spatial  and  Temporal  Evolution  of  the  Microwave  Streamer  at  (a)  0.0  ns, 
(b)  12.0  ns,  (c)  14.0  ns,  and  (d)  15.8  ns,  Respectively.  In  Each  Set  of  Plot,  the  Upper  Two 
Figures  Show  the  Electron  Density  and  the  Effective  Field  Distributions  in  a  2D  Cut,  and 
the  Lower  Two  Figures  Show  the  Corresponding  Distributions  as  ID  Curves  Along  the 

Center  Line  (the  Symmetric  Axis) 
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The  electron  density  and  the  secondary  electric  field  distribution  at  15.8  ns  are  shown  in  Figure 
16(a)  and  (b),  respectively.  From  Figure  16(a),  it  is  clear  that  the  plasma  is  elongated  in  the 
direction  of  the  incident  electric  field.  From  Figure  16(b),  the  plasma  shielding  to  the 
electromagnetic  waves  can  be  seen  clearly.  It  can  also  be  observed  that  the  secondary  electric 
field  is  concentrated  at  the  two  tips  of  the  plasma  bulk,  which  is  because  the  electric  field 
radiated  by  a  small  electric  dipole  (a  highly  oscillating  electron)  in  the  plasma  front  cannot  be 
cancelled  by  the  electric  fields  radiated  by  its  neighboring  dipoles  due  to  the  imbalanced  electron 
density  resulted  from  the  high  electron  density  gradient  at  the  plasma  front  [refer  to  Figure 
15(d)].  The  radiation  and  shielding  of  the  electromagnetic  fields  by  the  oscillating  plasma  can  be 
observed  in  Figure  16(c),  where  the  two  plots  are  presented  on  top  of  each  other. 


(c) 

Figure  16.  (a)  Electron  Density,  (b)  Secondary  Electric  Field  Distribution,  and  (c) 
Combination  View  of  the  Electron  Density  and  the  Secondary  Electric  Field  Generated  by 

the  Electron  Plasma  at  15.8  ns 


As  pointed  out  in  the  preceding  sections,  the  field  continuity  in  all  directions  is  very  important  in 
order  to  obtain  the  correct  plasma  coefficients  in  the  simulation.  This  can  be  guaranteed  by  the 
use  of  the  nodal  DGTD  method.  Presented  in  Figure  17  are  the  y  and  z  components  of  the 
secondary  electric  field  radiated  by  the  electron  plasma  at  15.8  ns,  which  demonstrate  the  field 
continuity  in  all  directions  throughout  the  simulation  domain. 
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(a)  (b) 

Figure  17.  The  (a)  y  and  (b)  z  Components  of  the  Secondary  Electric  Field  Radiated  by  the 
Electron  Plasma  at  15.8  ns,  Which  Demonstrate  the  Field  Continuity  in  all  Directions 
Throughout  the  Simulation  Domain.  The  x  Component  is  Zero  in  the  Plane  of  Observation 

and  is  Not  Shown 


Figure  18  demonstrates  the  change  of  polynomial  orders  at  the  initial  and  final  time  steps. 
Apparently,  the  proposed  dynamic  p-adaption  algorithm  can  track  the  plasma  front  when  it 
propagates  and  determine  the  appropriate  polynomial  order  for  different  positions  at  different 
time  steps.  With  the  aid  of  the  dynamic  /^-adaption  algorithm,  the  computational  time  can  be 
reduced  by  two  orders  of  magnitude  compared  with  that  using  a  uniform  polynomial  order 
throughout  the  entire  computational  domain. 


(a)  (b) 

Figure  18.  Dynamic  /^-Adaption  During  Simulation.  The  Upper  Half  of  the  Two  Figures 
Show  the  Polynomial  Orders  in  the  Computational  Domain,  the  Lower  Left  and  Lower 
Right  of  the  Two  Figures  are  the  Electron  Density  and  Electric  Field  Magnitude, 

Respectively,  (a)  0.0  ns;  (b)  15.8  ns 


5.0  CONCLUSIONS 

In  this  report,  a  highly  accurate  and  efficient  simulation  tool  is  developed  for  the  multi-physics 
and  multiscale  modeling  and  simulation  of  the  electromagnetic-plasma  interaction  and  the  high- 
power  microwave  breakdown  in  air.  The  coupled  electromagnetic-plasma  system  consists  of  two 
Maxwell’s  equations,  a  particle  kinetic  equation,  and  a  particle  diffusion  equation.  The  nonlinear 
dependence  of  the  plasma  parameters  on  the  electromagnetic  fields  makes  the  entire  system 
highly  nonlinear.  In  order  to  deal  with  such  a  multi-physics,  multiscale,  and  nonlinear  system, 
the  nodal  DGTD  method  is  adopted  together  with  the  high-order  Runge-Kutta  method  in  the 
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solution  of  the  hyperbolic  Maxwell’s  equations  and  the  parabolic  diffusion  equation.  To  obtain  a 
better  computational  efficiency,  a  dynamic  /^-adaption  algorithm  is  proposed  and  a  GPU- 
accelerated  algorithm  is  implemented.  With  the  development  of  all  the  advanced  computational 
techniques,  the  high-power  microwave  breakdown  problems  can  be  successfully  simulated, 
which  provides  deeper  insight  and  better  understanding  of  the  physical  phenomenon  and  the 
breakdown  process.  Further  effort  is  being  taken  to  test  on  the  nodal  DGTD  for  modeling  of 
various  devices  involving  dielectrics  and  conductors  where  a  normal  discontinuity  at  a  material 
interface  is  present. 
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